home *** CD-ROM | disk | FTP | other *** search
/ CD Ware Multimedia 1995 May / cd Ware (Juegos) Epimundo.iso / DOS / C / CMATH.ZIP / POW.C < prev    next >
Encoding:
C/C++ Source or Header  |  1989-03-21  |  11.4 KB  |  518 lines
  1. /*                            pow.c
  2.  *
  3.  *    Power function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, z, pow();
  10.  *
  11.  * z = pow( x, y );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Computes x raised to the yth power.  Analytically,
  18.  *
  19.  *      x**y  =  exp( y log(x) ).
  20.  *
  21.  * Following Cody and Waite, this program uses a lookup table
  22.  * of 2**-i/16 and pseudo extended precision arithmetic to
  23.  * obtain an extra three bits of accuracy in both the logarithm
  24.  * and the exponential.
  25.  *
  26.  *
  27.  *
  28.  * ACCURACY:
  29.  *
  30.  *                      Relative error:
  31.  * arithmetic   domain     # trials      peak         rms
  32.  *    IEEE     -26,26       30000      4.2e-16      7.7e-17
  33.  *    DEC      -26,26       18000      4.4e-17      8.9e-18
  34.  * 1/26 < x < 26, with log(x) uniformly distributed.
  35.  * -26 < y < 26, y uniformly distributed.
  36.  *    IEEE     0,8700       30000      1.5e-14      2.1e-15
  37.  *    DEC      0,8700        1000      1.3e-15      2.1e-16
  38.  * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
  39.  *
  40.  *
  41.  * ERROR MESSAGES:
  42.  *
  43.  *   message         condition      value returned
  44.  * pow overflow     x**y > MAXNUM      MAXNUM
  45.  * pow underflow   x**y < 1/MAXNUM       0.0
  46.  * pow domain      x<0 and y noninteger  0.0
  47.  *
  48.  */
  49.  
  50. /*
  51. Cephes Math Library Release 2.1:  December, 1988
  52. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  53. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  54. */
  55.  
  56.  
  57. #include "mconf.h"
  58. static char fname[] = {"pow"};
  59.  
  60. #define SQRTH 0.70710678118654752440
  61.  
  62. #ifdef UNK
  63. static double P[] = {
  64.   4.97778295871696322025E-1,
  65.   3.73336776063286838734E0,
  66.   7.69994162726912503298E0,
  67.   4.66651806774358464979E0
  68. };
  69. static double Q[] = {
  70. /* 1.00000000000000000000E0, */
  71.   9.33340916416696166113E0,
  72.   2.79999886606328401649E1,
  73.   3.35994905342304405431E1,
  74.   1.39995542032307539578E1
  75. };
  76. /* 2^(-i/16), IEEE precision */
  77. static double A[] = {
  78.   1.00000000000000000000E0,
  79.   9.57603280698573700036E-1,
  80.   9.17004043204671215328E-1,
  81.   8.78126080186649726755E-1,
  82.   8.40896415253714502036E-1,
  83.   8.05245165974627141736E-1,
  84.   7.71105412703970372057E-1,
  85.   7.38413072969749673113E-1,
  86.   7.07106781186547572737E-1,
  87.   6.77127773468446325644E-1,
  88.   6.48419777325504820276E-1,
  89.   6.20928906036742001007E-1,
  90.   5.94603557501360513449E-1,
  91.   5.69394317378345782288E-1,
  92.   5.45253866332628844837E-1,
  93.   5.22136891213706877402E-1,
  94.   5.00000000000000000000E-1
  95. };
  96. static double B[] = {
  97.  0.00000000000000000000E0,
  98.  1.64155361212281360176E-17,
  99.  4.09950501029074826006E-17,
  100.  3.97491740484881042808E-17,
  101. -4.83364665672645672553E-17,
  102.  1.26912513974441574796E-17,
  103.  1.99100761573282305549E-17,
  104. -1.52339103990623557348E-17,
  105.  0.00000000000000000000E0
  106. };
  107. static double R[] = {
  108.  1.49664108433729301083E-5,
  109.  1.54010762792771901396E-4,
  110.  1.33335476964097721140E-3,
  111.  9.61812908476554225149E-3,
  112.  5.55041086645832347466E-2,
  113.  2.40226506959099779976E-1,
  114.  6.93147180559945308821E-1
  115. };
  116.  
  117. #define douba(k) A[k]
  118. #define doubb(k) B[k]
  119. #define MEXP 2031.0
  120. #define MNEXP -2031.0
  121. #endif
  122.  
  123. #ifdef DEC
  124. static short P[] = {
  125. 0037776,0156313,0175332,0163602,
  126. 0040556,0167577,0052366,0174245,
  127. 0040766,0062753,0175707,0055564,
  128. 0040625,0052035,0131344,0155636,
  129. };
  130. static short Q[] = {
  131. /*0040200,0000000,0000000,0000000,*/
  132. 0041025,0052644,0154404,0105155,
  133. 0041337,0177772,0007016,0047646,
  134. 0041406,0062740,0154273,0020020,
  135. 0041137,0177054,0106127,0044555,
  136. };
  137. static short A[] = {
  138. 0040200,0000000,0000000,0000000,
  139. 0040165,0022575,0012444,0103314,
  140. 0040152,0140306,0163735,0022071,
  141. 0040140,0146336,0166052,0112341,
  142. 0040127,0042374,0145326,0116553,
  143. 0040116,0022214,0012437,0102201,
  144. 0040105,0063452,0010525,0003333,
  145. 0040075,0004243,0117530,0006067,
  146. 0040065,0002363,0031771,0157145,
  147. 0040055,0054076,0165102,0120513,
  148. 0040045,0177326,0124661,0050471,
  149. 0040036,0172462,0060221,0120422,
  150. 0040030,0033760,0050615,0134251,
  151. 0040021,0141723,0071653,0010703,
  152. 0040013,0112701,0161752,0105727,
  153. 0040005,0125303,0063714,0044173,
  154. 0040000,0000000,0000000,0000000
  155. };
  156. static short B[] = {
  157. 0000000,0000000,0000000,0000000,
  158. 0021473,0040265,0153315,0140671,
  159. 0121074,0062627,0042146,0176454,
  160. 0121413,0003524,0136332,0066212,
  161. 0121767,0046404,0166231,0012553,
  162. 0121257,0015024,0002357,0043574,
  163. 0021736,0106532,0043060,0056206,
  164. 0121310,0020334,0165705,0035326,
  165. 0000000,0000000,0000000,0000000
  166. };
  167.  
  168. static short R[] = {
  169. 0034173,0014076,0137624,0115771,
  170. 0035041,0076763,0003744,0111311,
  171. 0035656,0141766,0041127,0074351,
  172. 0036435,0112533,0073611,0116664,
  173. 0037143,0054106,0134040,0152223,
  174. 0037565,0176757,0176026,0025551,
  175. 0040061,0071027,0173721,0147572
  176. };
  177.  
  178. /*
  179. static double R[] = {
  180. 0.14928852680595608186e-4,
  181. 0.15400290440989764601e-3,
  182. 0.13333541313585784703e-2,
  183. 0.96181290595172416964e-2,
  184. 0.55504108664085595326e-1,
  185. 0.24022650695909537056e0,
  186. 0.69314718055994529629e0
  187. };
  188. */
  189. #define douba(k) (*(double *)&A[(k)<<2])
  190. #define doubb(k) (*(double *)&B[(k)<<2])
  191. #define MEXP 2031.0
  192. #define MNEXP -2031.0
  193. #endif
  194.  
  195. #ifdef IBMPC
  196. static short P[] = {
  197. 0x5cf0,0x7f5b,0xdb99,0x3fdf,
  198. 0xdf15,0xea9e,0xddef,0x400d,
  199. 0xeb6f,0x7f78,0xccbd,0x401e,
  200. 0x9b74,0xb65c,0xaa83,0x4012,
  201. };
  202. static short Q[] = {
  203. /*0x0000,0x0000,0x0000,0x3ff0,*/
  204. 0x914e,0x9b20,0xaab4,0x4022,
  205. 0xc9f5,0x41c1,0xffff,0x403b,
  206. 0x6402,0x1b17,0xccbc,0x4040,
  207. 0xe92e,0x918a,0xffc5,0x402b,
  208. };
  209. static short A[] = {
  210. 0x0000,0x0000,0x0000,0x3ff0,
  211. 0x90da,0xa2a4,0xa4af,0x3fee,
  212. 0xa487,0xdcfb,0x5818,0x3fed,
  213. 0x529c,0xdd85,0x199b,0x3fec,
  214. 0xd3ad,0x995a,0xe89f,0x3fea,
  215. 0xf090,0x82a3,0xc491,0x3fe9,
  216. 0xa0db,0x422a,0xace5,0x3fe8,
  217. 0x0187,0x73eb,0xa114,0x3fe7,
  218. 0x3bcd,0x667f,0xa09e,0x3fe6,
  219. 0x5429,0xdd48,0xab07,0x3fe5,
  220. 0x2a27,0xd536,0xbfda,0x3fe4,
  221. 0x3422,0x4c12,0xdea6,0x3fe3,
  222. 0xb715,0x0a31,0x06fe,0x3fe3,
  223. 0x6238,0x6e75,0x387a,0x3fe2,
  224. 0x517b,0x3c7d,0x72b8,0x3fe1,
  225. 0x890f,0x6cf9,0xb558,0x3fe0,
  226. 0x0000,0x0000,0x0000,0x3fe0
  227. };
  228. static short B[] = {
  229. 0x0000,0x0000,0x0000,0x0000,
  230. 0x3707,0xd75b,0xed02,0x3c72,
  231. 0xcc81,0x345d,0xa1cd,0x3c87,
  232. 0x4b27,0x5686,0xe9f1,0x3c86,
  233. 0x6456,0x13b2,0xdd34,0xbc8b,
  234. 0x42e2,0xafec,0x4397,0x3c6d,
  235. 0x82e4,0xd231,0xf46a,0x3c76,
  236. 0x8a76,0xb9d7,0x9041,0xbc71,
  237. 0x0000,0x0000,0x0000,0x0000
  238. };
  239. static short R[] = {
  240. 0x937f,0xd7f2,0x6307,0x3eef,
  241. 0x9259,0x60fc,0x2fbe,0x3f24,
  242. 0xef1d,0xc84a,0xd87e,0x3f55,
  243. 0x33b7,0x6ef1,0xb2ab,0x3f83,
  244. 0x1a92,0xd704,0x6b08,0x3fac,
  245. 0xc56d,0xff82,0xbfbd,0x3fce,
  246. 0x39ef,0xfefa,0x2e42,0x3fe6
  247. };
  248.  
  249. #define douba(k) (*(double *)&A[(k)<<2])
  250. #define doubb(k) (*(double *)&B[(k)<<2])
  251. #define MEXP 16383.0
  252. /* The following if denormal numbers are supported, else -MEXP: */
  253. #define MNEXP -17183.0
  254. #endif
  255.  
  256. #ifdef MIEEE
  257. static short P[] = {
  258. 0x3fdf,0xdb99,0x7f5b,0x5cf0,
  259. 0x400d,0xddef,0xea9e,0xdf15,
  260. 0x401e,0xccbd,0x7f78,0xeb6f,
  261. 0x4012,0xaa83,0xb65c,0x9b74
  262. };
  263. static short Q[] = {
  264. 0x4022,0xaab4,0x9b20,0x914e,
  265. 0x403b,0xffff,0x41c1,0xc9f5,
  266. 0x4040,0xccbc,0x1b17,0x6402,
  267. 0x402b,0xffc5,0x918a,0xe92e
  268. };
  269. static short A[] = {
  270. 0x3ff0,0x0000,0x0000,0x0000,
  271. 0x3fee,0xa4af,0xa2a4,0x90da,
  272. 0x3fed,0x5818,0xdcfb,0xa487,
  273. 0x3fec,0x199b,0xdd85,0x529c,
  274. 0x3fea,0xe89f,0x995a,0xd3ad,
  275. 0x3fe9,0xc491,0x82a3,0xf090,
  276. 0x3fe8,0xace5,0x422a,0xa0db,
  277. 0x3fe7,0xa114,0x73eb,0x0187,
  278. 0x3fe6,0xa09e,0x667f,0x3bcd,
  279. 0x3fe5,0xab07,0xdd48,0x5429,
  280. 0x3fe4,0xbfda,0xd536,0x2a27,
  281. 0x3fe3,0xdea6,0x4c12,0x3422,
  282. 0x3fe3,0x06fe,0x0a31,0xb715,
  283. 0x3fe2,0x387a,0x6e75,0x6238,
  284. 0x3fe1,0x72b8,0x3c7d,0x517b,
  285. 0x3fe0,0xb558,0x6cf9,0x890f,
  286. 0x3fe0,0x0000,0x0000,0x0000
  287. };
  288. static short B[] = {
  289. 0x0000,0x0000,0x0000,0x0000,
  290. 0x3c72,0xed02,0xd75b,0x3707,
  291. 0x3c87,0xa1cd,0x345d,0xcc81,
  292. 0x3c86,0xe9f1,0x5686,0x4b27,
  293. 0xbc8b,0xdd34,0x13b2,0x6456,
  294. 0x3c6d,0x4397,0xafec,0x42e2,
  295. 0x3c76,0xf46a,0xd231,0x82e4,
  296. 0xbc71,0x9041,0xb9d7,0x8a76,
  297. 0x0000,0x0000,0x0000,0x0000
  298. };
  299. static short R[] = {
  300. 0x3eef,0x6307,0xd7f2,0x937f,
  301. 0x3f24,0x2fbe,0x60fc,0x9259,
  302. 0x3f55,0xd87e,0xc84a,0xef1d,
  303. 0x3f83,0xb2ab,0x6ef1,0x33b7,
  304. 0x3fac,0x6b08,0xd704,0x1a92,
  305. 0x3fce,0xbfbd,0xff82,0xc56d,
  306. 0x3fe6,0x2e42,0xfefa,0x39ef
  307. };
  308.  
  309. #define douba(k) (*(double *)&A[(k)<<2])
  310. #define doubb(k) (*(double *)&B[(k)<<2])
  311. #define MEXP 16383.0
  312. #define MNEXP -17183.0
  313. #endif
  314.  
  315. /* log2(e) - 1 */
  316. #define LOG2EA 0.44269504088896340736
  317.  
  318. #define F W
  319. #define Fa Wa
  320. #define Fb Wb
  321. #define G W
  322. #define Ga Wa
  323. #define Gb u
  324. #define H W
  325. #define Ha Wb
  326. #define Hb Wb
  327.  
  328. extern double MAXNUM;
  329.  
  330. double pow( x, y )
  331. double x, y;
  332. {
  333. double w, z, W, Wa, Wb, ya, yb;
  334. double u, v;
  335. /* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
  336. int e, i, nflg;
  337. double floor(), fabs(), frexp(), ldexp();
  338. double reduc(), polevl(), p1evl(), powi();
  339.  
  340.  
  341. nflg = 0;    /* flag = 1 if x<0 raised to integer power */
  342. w = floor(y);
  343. if( (w == y) && (fabs(w) < 32768.0) )
  344.     {
  345.     i = w;
  346.     w = powi( x, i );
  347.     return( w );
  348.     }
  349.  
  350.  
  351. if( x <= 0.0 )
  352.     {
  353.     if( x == 0.0 )
  354.         {
  355.         if( y == 0.0 )
  356.             return( 1.0 );  /*   0**0   */
  357.         else  
  358.             return( 0.0 );  /*   0**y   */
  359.         }
  360.     else
  361.         {
  362.         if( w != y )
  363.             { /* noninteger power of negative number */
  364.             mtherr( fname, DOMAIN );
  365.             return(0.0);
  366.             }
  367.         nflg = 1;
  368.         x = fabs(x);
  369.         }
  370.     }
  371.  
  372. /* separate significand from exponent */
  373. x = frexp( x, &e );
  374.  
  375. /* find significand in antilog table A[] */
  376. i = 1;
  377. if( x <= douba(9) )
  378.     i = 9;
  379. if( x <= douba(i+4) )
  380.     i += 4;
  381. if( x <= douba(i+2) )
  382.     i += 2;
  383. if( x >= douba(1) )
  384.     i = -1;
  385. i += 1;
  386.  
  387.  
  388. /* Find (x - A[i])/A[i]
  389.  * in order to compute log(x/A[i]):
  390.  *
  391.  * log(x) = log( a x/a ) = log(a) + log(x/a)
  392.  *
  393.  * log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a
  394.  */
  395. x -= douba(i);
  396. x -= doubb(i/2);
  397. x /= douba(i);
  398.  
  399.  
  400. /* rational approximation for log(1+v):
  401.  *
  402.  * log(1+v)  =  v  -  v**2/2  +  v**3 P(v) / Q(v)
  403.  */
  404. z = x*x;
  405. w = x * ( z * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) );
  406. w = w - ldexp( z, -1 );   /*  w - 0.5 * z  */
  407.  
  408. /* Convert to base 2 logarithm:
  409.  * multiply by log2(e)
  410.  */
  411. w = w + LOG2EA * w;
  412. /* Note x was not yet added in
  413.  * to above rational approximation,
  414.  * so do it now, while multiplying
  415.  * by log2(e).
  416.  */
  417. z = w + LOG2EA * x;
  418. z = z + x;
  419.  
  420. /* Compute exponent term of the base 2 logarithm. */
  421. w = -i;
  422. w = ldexp( w, -4 );    /* divide by 16 */
  423. w += e;
  424. /* Now base 2 log of x is w + z. */
  425.  
  426. /* Multiply base 2 log by y, in extended precision.
  427.  
  428. /* separate y into large part ya
  429.  * and small part yb less than 1/16
  430.  */
  431. ya = reduc(y);
  432. yb = y - ya;
  433.  
  434.  
  435. F = z * y  +  w * yb;
  436. Fa = reduc(F);
  437. Fb = F - Fa;
  438.  
  439. G = Fa + w * ya;
  440. Ga = reduc(G);
  441. Gb = G - Ga;
  442.  
  443. H = Fb + Gb;
  444. Ha = reduc(H);
  445. w = ldexp( Ga+Ha, 4 );
  446.  
  447. /* Test the power of 2 for overflow */
  448. if( w > MEXP )
  449.     {
  450.     mtherr( fname, OVERFLOW );
  451.     return( MAXNUM );
  452.     }
  453.  
  454. if( w < MNEXP )
  455.     {
  456.     mtherr( fname, UNDERFLOW );
  457.     return( 0.0 );
  458.     }
  459.  
  460. e = w;
  461. Hb = H - Ha;
  462.  
  463. if( Hb > 0.0 )
  464.     {
  465.     e += 1;
  466.     Hb -= 0.0625;
  467.     }
  468.  
  469. /* Now the product y * log2(x)  =  Hb + e/16.0.
  470.  *
  471.  * Compute base 2 exponential of Hb,
  472.  * where -0.0625 <= Hb <= 0.
  473.  */
  474. z = Hb * polevl( Hb, R, 6 );  /*    z  =  2**Hb - 1    */
  475.  
  476. /* Express e/16 as an integer plus a negative number of 16ths.
  477.  * Find lookup table entry for the fractional power of 2.
  478.  */
  479. if( e < 0 )
  480.     i = 0;
  481. else
  482.     i = 1;
  483. i = e/16 + i;
  484. e = 16*i - e;
  485. w = douba( e );
  486. z = w + w * z;      /*    2**-e * ( 1 + (2**Hb-1) )    */
  487. z = ldexp( z, i );  /* multiply by integer power of 2 */
  488.  
  489. if( nflg )
  490.     {
  491. /* For negative x,
  492.  * find out if the integer exponent
  493.  * is odd or even.
  494.  */
  495.     w = ldexp( y, -1 );
  496.     w = floor(w);
  497.     w = ldexp( w, 1 );
  498.     if( w != y )
  499.         z = -z; /* odd exponent */
  500.     }
  501.  
  502. return( z );
  503. }
  504.  
  505.  
  506. /* Find a multiple of 1/16 that is within 1/16 of x. */
  507. static double reduc(x)
  508. double x;
  509. {
  510. double t;
  511. double ldexp(), floor();
  512.  
  513. t = ldexp( x, 4 );
  514. t = floor( t );
  515. t = ldexp( t, -4 );
  516. return(t);
  517. }
  518.